*************************************************************************************;
*************************************************************************************;
*************************************************************************************;
***                                                                               ***;
*** ICCs Macro (SAS Version 9.1) -- Most recent update, July 25, 2005             ***;
*** Written by -- Paul Weiss and George Cotsonis                                  ***;
*** Under the Direction of -- Michael Haber                                       ***;
*** Supported by NIH award R01 MH070028                                           ***;
***                                                                               ***;
*** This macro will generate measures of agreement and consistency between raters ***;
*** using moments and components of variance. The macro requires that the data be ***;
*** in columnar form, with each individual observation representing a case, and   ***;
*** each column representing the measurements taken by an individual rater.       ***;
*** Each case needs a unique identifier, called ID, which must appear in the      ***;
*** data set as a variable.                                                       ***;
***                                                                               ***;
***                                                                               ***;
*** Example:                                                                      ***;
***                                                                               ***;
*** ID Measure1 Measure2 Measure3 Measure4                                        ***;
***  1    8         8        9        9                                           ***;
***  2    9         9        8        8                                           ***;
***  3    8         9        9        8                                           ***;
***  4    9         8        8        9                                           ***;
***                                                                               ***;
*** The macro requires the user to specify the following parameters:              ***;
***                                                                               ***;
*** lib     : the library where the data are stored (WORK is default)             ***;
*** dsn     : the data set for analysis                                           ***;
*** rep     : replicate variable name [Not Working]                               ***;
*** var     : a one-word description of the measurement                           ***;
*** nummeas : the number of measurers                                             ***;
*** numsub  : the number of individuals being measured                            ***;
*** numrep  : the number of replicates on each individual (default is 1)          ***;
*** varlist : an explicit listing of variables for the analysis                   ***;
*** dataname: explicit description of data source                                 ***;
*** seed    : random number seed for method 2 (random selection of measure)       ***;
*** alpha   : significance level for confidence intervals                         ***;
*** idvar   : the variable indicating the unique ID for each subject              ***;
*** boot    : the number of bootstrap samples desired (0 suppresses this feature) ***;
***                                                                               ***;
*************************************************************************************;
***                                                                               ***;
*** Due to confusion on varibles for subject identifiers, the new version allows  ***;
*** the user to define the subject identifier, and this is recoded to "id" in     ***;
*** the macro to prevent errors from naming issues arising from capitalization.   ***;
***                                                                               ***;
*************************************************************************************;
*************************************************************************************;

%macro ICCs (lib=, dsn=, rep=, var=, numrater=, nummeas=, numsub=, numrep=, varlist=, dataname=, 
             seed=0, alpha=0.05, idvar=id, boot=0 ) ;



%let alphalow = %sysevalf(0.5*&alpha);
%let alphahi = %sysevalf(1-(0.5*&alpha));
%let conf = %sysevalf(100*(1-&alpha));



footnote;

*** Page 1: First 25 observations in the data ***;

Title1 ;
title2 ;
title3 "&dataname -- &var";

%if &numsub gt 25 %then %do;
Title4 "First 25 Observations in the Data";
%end; %else %do;
Title4 "First &numsub Observations in the Data";
%end;

ods select all;

proc print data=&lib..&dsn (obs=25);
   var &idvar &varlist;
run;

%if &numrep eq 1 %then %do; *** Begin NON-REPLICATE ICC calculations ***;

*** Page 2: Means And Standard Deviations for each observer ***;

Title1 ;
title2 ;
title3 "&dataname -- &var";
Title4 "Simple Descriptive Statistics For Each Observer";

proc means mean std data=&lib..&dsn noprint;
   var &varlist;
   output out=means mean=
   %do bar=1 %to &nummeas;
   mean&var&bar
   %end;
   std=
   %do sd=1 %to &nummeas;
   sd&var&sd
   %end;
   ;
run;

%do ds=1 %to &nummeas;

data rater&ds;
   set means;
   mean=mean&var&ds;
   sd  =sd&var&ds;
   keep mean sd;
run;

%end;

data raters;
   set
   %do ds=1 %to &nummeas;
   rater&ds
   %end;
   ;
   rater+1;
   label rater="Observer"
         mean="Mean"
         sd="Standard Deviation"
   ;
run;

proc print noobs label;
   var rater mean sd;
run;


*** Start with Pearson Correlations ***;

Title1 ;
title2 ;
title3 "&dataname -- &var";
Title4 "Pearson Correlation Coefficients and Covariances";

*** Preparing data for analysis ***;

ods select none;
proc corr data=&lib..&dsn cov;
   var &varlist;
   ods output cov=covmat;
run;

%do d=1 %to &nummeas;

data a&d;
   set covmat;
   if _n_=&d;
   %do j=1 %to &nummeas;
   &var&d&j=&var&j;
   %end;
run;

%end;

data oneline;
   merge
   %do d=1 %to &nummeas;
   a&d
   %end;
   ;;
   keep
   %do i=1 %to &nummeas;
       %do j=1 %to &nummeas;
	       &var&i&j
       %end;
   %end;
   ;; 
run;

data pear;
   set oneline;

   %do i=1 %to &nummeas;
       %do j=1 %to &nummeas;

	   r&i&j = &var&i&j / sqrt(&var&i&i * &var&j&j);

	   %end;
   %end;

run;

%let a11 = 11;

data corrs;
   set pear;
   array r(&nummeas,&nummeas) r11 -- r&nummeas&nummeas;
   array x(&nummeas,&nummeas) &var&a11 -- &var&nummeas&nummeas;
   do i=1 to &nummeas;
      do j=1 to &nummeas;
	     pearson=r(i,j);
		 cov=x(i,j);
		 var1=x(i,i);
		 var2=x(j,j);
		 if i lt j then output;
	  end;
   end;
   keep i j pearson cov var1 var2;
   label pearson="Pearson Correlation"
         cov="Covariance"
         var1="Variance 1"
		 var2="Variance 2"
		 i="Measure 1"
		 j="Measure 2"
		 ;
run;

ods select all;
proc print label noobs;
   var i j pearson cov ;
run;

*** Calculating the OCCC ***;

Title4 "Coefficients of Agreement";

* Need Covariances from Pearson *;

proc means sum noprint data=corrs;
   var cov;
   output out=covs sum=sumcov;
run;

proc sort data=corrs;
   by i j;
run;

data v1;
   set corrs;
   by i j;
   if first.i; *** Keeps first j-1 variances ***;
   measure=i;
   var=var1;
   keep measure var;
run;
   
data v2;
   set corrs;
   by i j;
   if last.j; *** Keeps last j-1 variances ***;
   measure=j;
   var=var2;
   keep measure var;
run;

data v;
   set v1 v2;
run;

proc sort;
   by measure;
run;

data vx;
   set v;
   by measure;
   if first.measure;
run;

proc means sum noprint data=vx;
   var var;
   output out=vars sum=sumvar;
run;

proc means mean data=&lib..&dsn noprint;
   var &varlist;
   output out=means mean=;
run;

data meandiffs;
   set means;
   array means (&nummeas) &varlist;
   xend = &nummeas - 1;
   meansq=0; * to start *;
   do i = 1 to xend;
      do j = i+1 to &nummeas;
	     meansq + (means(i) - means(j))**2;
	  end;
   end;
   drop i j xend;
run;

proc means mean data=corrs noprint;
   var pearson;
   output out=meancorr mean=corrbar;
run;

data occc;
   merge covs vars meandiffs meancorr;
   nab = 2*sumcov / (((&nummeas - 1)*sumvar) + meansq);
   nad = 2*sumcov / (((&nummeas - 1)*sumvar));
   nl  = corrbar;
   label nab = "Coefficient of Absolute Agreement (CCC)"
         nad = "Coefficient of Additive Agreement"
		 nl  = "Coefficient of Linear Agreement"
         sumcov = "Sum of Covariance"
		 sumvar = "Sum of Variance"
		 meansq  = "Sum of Squared Difference in Means"
		 ;
proc print label noobs;
   var sumcov sumvar meansq nab nad nl;
run;

*** pairwise correlations for absolute and additive agreement ***;

%if &nummeas ge 3 %then %do; 

%do x=1 %to %eval(&nummeas - 1);
%do y=%eval(&x+1) %to &nummeas;

data posthoc; *** Create a data set with only two measures, repeat for all pairwise measures ***;
   set &lib..&dsn;
   v1=&var&x; * Pick first measure *;
   v2=&var&y; * Pick second measure *;
   keep id v1 v2; * Keep data set with chosen measures as measure 1 and measure 2 *;
run;

ods select none;
proc corr data=posthoc cov;
   var v1 v2;
   ods output cov=covmat;
run;

%do d=1 %to 2;

data a&d;
   set covmat;
   if _n_=&d;
   %do j=1 %to 2;
   &var&d&j=v&j;
   %end;
run;

%end;

data oneline;
   merge
   %do d=1 %to 2;
   a&d
   %end;
   ;;
   keep
   %do i=1 %to 2;
       %do j=1 %to 2;
	       &var&i&j
       %end;
   %end;
   ;; 
run;

data pear;
   set oneline;

   %do i=1 %to 2;
       %do j=1 %to 2;

	   r&i&j = &var&i&j / sqrt(&var&i&i * &var&j&j);

	   %end;
   %end;

run;

%let a11 = 11;
%let a22 = 22;

data corrs;
   set pear;
   array r(2,2) r11 -- r22;
   array x(2,2) &var&a11 -- &var&a22;
   do i=1 to 2;
      do j=1 to 2;
	     pearson=r(i,j);
		 cov=x(i,j);
		 var1=x(i,i);
		 var2=x(j,j);
		 if i lt j then output;
	  end;
   end;
   keep i j pearson cov var1 var2;
   label pearson="Pearson Correlation"
         cov="Covariance"
         var1="Variance 1"
		 var2="Variance 2"
		 i="Measure 1"
		 j="Measure 2"
		 ;
run;


proc means mean data=posthoc noprint;
   var v1 v2;
   output out=means mean=vb1 vb2;
run;

data meandiffs;
   set means;
   meansq = (vb1 - vb2)**2;
run;

proc means mean data=corrs noprint;
   var pearson;
   output out=meancorr mean=corrbar;
run;

data occc&x&y;
   merge corrs meancorr meandiffs ;
   rater1=&x;
   rater2=&y;
   nab = 2*cov / (var1 + var2 + meansq);
   nad = 2*cov / (var1 + var2);
   nl  = corrbar;
   label nab = "Coefficient of Absolute Agreement (CCC)"
         nad = "Coefficient of Additive Agreement"
		 nl  = "Coefficient of Linear Agreement"
         cov = "Covariance"
		 var = "Variance"
		 meansq  = "Squared Difference in Means"
		 ;
run;

%end;
%end;

data occcpost;
   set
     %do x=1 %to %eval(&nummeas - 1);
        %do y=%eval(&x+1) %to &nummeas;
           occc&x&y
	    %end;
     %end;
   ;;
run;

ods select all;

proc print label noobs;
   var rater1 rater2 nab nad nl;
   title5 "Coefficients of Agreement, Pairwise by Rater";
run;

ods select none;

%end;



*** Calculating ICC-C and ICC-A ***;

Title4 "Intraclass Correlation Coefficients for Agreement and Consistency";


data icc;

   set &lib..&dsn;

   %do i=1 %to &nummeas;

       id=&idvar; &var=&var&i; rater=&i; output;

   %end;
   keep id &var rater;

run;

ods select none;
proc glm data=icc;
   class &idvar rater ;
   model &var = &idvar rater ;
   random &idvar ;
   ods output OverallANOVA=over ModelANOVA=SS ;
run;

data one;
   set over;
   if Source="Error";
   edf=df; ems=ms;
   keep edf ems;
run;

data two;
   set ss;
   if Hypothesistype=3;
   if source="&idvar" then do;
      bdf=df; bms=ms; jdf=0; jms=0; fobs=fvalue;
   end; else
   if source='rater' then do;
      bdf=0; bms=0; jdf=df; jms=ms; fobs=0;
   end;

proc means sum noprint data=two;
   var bdf bms jdf jms fobs;
   output out=twoa sum= ;
run;

data iccfull;
   merge one twoa;
   drop _type_ _freq_;

   n=&numrep;    * This will be assigned in a macro variable later as number of replicates *;
   a=&nummeas;   * This will be assigned in a macro variable later as number of raters     *;
   b=&numsub;    * This will be assigned in a macro variable later as number of subjects   *;

   sse = ems;

   sss = (bms - ems)/(n * a);

   ssr = (jms - ems)/(n * b);

   ICCa = sss / (sss + ssr + sse);
   ICCc = sss / (sss + sse);

   conf_lim=1-(&alpha/2);

   *** ICC-C Confidence Limits ***;
   cdfnum=b-1;
   cdfden=(b-1)*(a-1);

   cftabledl=finv(conf_lim, cdfnum, cdfden);
   cftabledu=finv(conf_lim, cdfden, cdfnum);

   fl = fobs/cftabledl;
   fu = fobs*cftabledu;

   ccll= (fl-1)/(fl + a - 1);
   cclu= (fu-1)/(fu + a - 1);

   *** Confidence Limits for ICC-A ***;

   aa=(a*ICCa)/(b*(1-ICCa));
   ab=1+((a*ICCa*(b-1))/(b*(1-ICCa)));

   asatnum=(aa*jms+ab*ems)**2;
   asatden1=((aa*jms)**2)/(a-1);
   asatden2=((ab*ems)**2)/((a-1)*(b-1));

   anu = asatnum/(asatden1 + asatden2);

   aftabledl=finv(conf_lim, b-1, anu);
   aftabledu=finv(conf_lim, anu, b-1);

   aclltop=b*(bms-aftabledl*ems);
   acllbot=aftabledl*(a*jms+((a*b) - a - b)*ems)+b*bms;

   acll=aclltop/acllbot;

   aclutop=b*(aftabledu*bms-ems);
   aclubot=a*jms+((a*b) - a - b)*ems + b*aftabledu*bms;

   aclu=aclutop/aclubot;

   rater1="Overall";
   rater2=" ";

   label sse="Error Component"
         ssr="Rater Component"
		 sss="Subject Component"
		 ICCa="Intraclass Correlation (Agreement)"
		 ICCc="Intraclass Correlation (Consistency)"
		 ccll="Lower Confidence Limit for ICC-C (alpha=&alpha)"
		 cclu="Upper Confidence Limit for ICC-C (alpha=&alpha)"
		 acll="Lower Confidence Limit for ICC-A (alpha=&alpha)"
		 aclu="Upper Confidence Limit for ICC-A (alpha=&alpha)"
   ;

run;



%if &nummeas ge 3 %then %do; *** Begin Post-hocs block ***;

%do x=1 %to %eval(&nummeas - 1);
%do y=%eval(&x+1) %to &nummeas;

data posthoc; *** Create a data set with only two measures, repeat for all pairwise measures ***;
   set &lib..&dsn;
   var1=&var&x; * Pick first measure *;
   var2=&var&y; * Pick second measure *;
   keep id var1 var2; * Keep data set with chosen measures as measure 1 and measure 2 *;
run;


data icc;

   set posthoc;

   %do i=1 %to 2;

       id=&idvar; var=var&i; rater=&i; output;

   %end;
   keep id var rater;

run;

ods select none;
proc glm data=icc;
   class id rater ;
   model var = id rater ;
   random id ;
   ods output OverallANOVA=over ModelANOVA=SS ;
run;

data one;
   set over;
   if Source="Error";
   edf=df; ems=ms;
   keep edf ems;
run;

data two;
   set ss;
   if Hypothesistype=3;
   if source="&idvar" then do;
      bdf=df; bms=ms; jdf=0; jms=0; fobs=fvalue;
   end; else
   if source='rater' then do;
      bdf=0; bms=0; jdf=df; jms=ms; fobs=0;
   end;

proc means sum noprint data=two;
   var bdf bms jdf jms fobs;
   output out=twoa sum= ;
run;

data iccpost&x&y;
   merge one twoa;
   drop _type_ _freq_;

   n=&numrep;    * This will be assigned in a macro variable later as number of replicates *;
   a=2       ;   * Use two (2) for pairwise post-hoc analysis in ICC calculations          *;
   b=&numsub;    * This will be assigned in a macro variable later as number of subjects   *;

   sse = ems;

   sss = (bms - ems)/(n * a);

   ssr = (jms - ems)/(n * b);

   ICCa = sss / (sss + ssr + sse);
   ICCc = sss / (sss + sse);

   conf_lim=1-(&alpha/2);

   *** ICC-C Confidence Limits ***;
   cdfnum=b-1;
   cdfden=(b-1)*(a-1);

   cftabledl=finv(conf_lim, cdfnum, cdfden);
   cftabledu=finv(conf_lim, cdfden, cdfnum);

   fl = fobs/cftabledl;
   fu = fobs*cftabledu;

   ccll= (fl-1)/(fl + a - 1);
   cclu= (fu-1)/(fu + a - 1);

   *** Confidence Limits for ICC-A ***;

   aa=(a*ICCa)/(b*(1-ICCa));
   ab=1+((a*ICCa*(b-1))/(b*(1-ICCa)));

   asatnum=(aa*jms+ab*ems)**2;
   asatden1=((aa*jms)**2)/(a-1);
   asatden2=((ab*ems)**2)/((a-1)*(b-1));

   anu = asatnum/(asatden1 + asatden2);

   aftabledl=finv(conf_lim, a-1, anu);
   aftabledu=finv(conf_lim, anu, a-1);

   aclltop=b*(bms-aftabledl*ems);
   acllbot=aftabledl*(b*jms+(ab - a -b)*ems)+b*bms;

   acll=aclltop/acllbot;

   aclutop=b*(aftabledu*bms-ems);
   aclubot=b*jms+(ab - a -b)*ems + b*aftabledu*bms;

   aclu=aclutop/aclubot;

   rater1="&x";
   rater2="&y";

   label sse="Error Component"
         ssr="Rater Component"
		 sss="Subject Component"
		 ICCa="Intraclass Correlation (Agreement)"
		 ICCc="Intraclass Correlation (Consistency)"
		 ccll="Lower Confidence Limit for ICC-C (alpha=&alpha)"
		 cclu="Upper Confidence Limit for ICC-C (alpha=&alpha)"
		 acll="Lower Confidence Limit for ICC-A (alpha=&alpha)"
		 aclu="Upper Confidence Limit for ICC-A (alpha=&alpha)"
		 rater1="Rater 1"
		 rater2="Rater 2"
   ;

run;
%end;
%end;
%end; *** End post-hocs block ***;

data icc;
   set iccfull 
   %if &nummeas ge 3 %then %do;
      %do x=1 %to %eval(&nummeas - 1);
          %do y=%eval(&x+1) %to &nummeas;
   iccpost&x&y
          %end;
      %end;
   %end;
   ;;
run;

ods select all;
title5;

proc print label noobs data=icc;
   var %if &nummeas ge 3 %then %do;
       rater1 rater2 sss ssr sse ;
       %end; %else %do;
       sss ssr sse ;
	   %end;
	   ;;
Title4 "ANOVA Variance Components";
run;

proc print label noobs data=icc;
   var %if &nummeas ge 3 %then %do;
       rater1 rater2 icca acll aclu ;
       %end; %else %do;
       icca acll aclu ;
	   %end;
	   ;;
Title4 "Intraclass Correlation Coefficients for Agreement";
run;

proc print label noobs data=icc;
   var %if &nummeas ge 3 %then %do;
       rater1 rater2 iccc ccll cclu;
       %end; %else %do;
       iccc ccll cclu;
	   %end;
	   ;;
Title4 "Intraclass Correlation Coefficients for Consistency";
run;


*** Other Measures of Pairwise Observer Agreement ***;

%if &nummeas ge 3 %then %do; *** Begin pairwise block ***;

%do x=1 %to %eval(&nummeas - 1);
%do y=%eval(&x+1) %to &nummeas;

data pairwise; *** Create a data set with only two measures, repeat for all pairwise measures ***;
   set &lib..&dsn;
   var1=&var&x; * Pick first measure *;
   var2=&var&y; * Pick second measure *;

   diff=var1 - var2;
   dsq=diff*diff;
   dabs=abs(diff);

   keep id var1 var2 diff dsq dabs; * Keep data set with chosen measures as measure 1 and measure 2 *;
run;

proc univariate noprint;
   var dabs dsq diff;
   output out=stats mean=mad msd mdiff
                    pctlpre=p_ 
                    pctlpts=80 90 
                    pctlname=abtdi80 abtdi90
   ;
run;

data stats&x&y;
   set stats;
   rater1=&x;
   rater2=&y;
run;

%end;
%end;

data omoa;
   set

     %do x=1 %to %eval(&nummeas - 1);
         %do y=%eval(&x+1) %to &nummeas;
             stats&x&y
		 %end;
	 %end;

   ;;
   label msd="Mean Squared Deviation"
         mad="Mean Absolute Deviation"
		 mdiff="Mean Difference"
		 P_abtdi80="TDI 80"
		 P_abtdi90="TDI 90"
		 rater1="Rater 1"
		 rater2="Rater 2"

		 ;
run;

ods select all;

proc print label noobs;
   var rater1 rater2 mdiff msd mad P_abtdi80 P_abtdi90;
   Title4 "Other Measures of Pairwise Observer Agreement";
run;

ods select none;

%end;  %else %do; *** Begin Pairwise block ***;

data pairwise; *** Create a data set with only two measures, repeat for all pairwise measures ***;
   set &lib..&dsn;
   array rater (2) &varlist;
   diff=rater(1) - rater(2);
   dsq=diff*diff;
   dabs=abs(diff);

   keep id &varlist diff dsq dabs; * Keep data set with chosen measures as measure 1 and measure 2 *;
run;

proc univariate noprint;
   var dabs dsq diff;
   output out=stats mean=mad msd mdiff
                    pctlpre=p_ 
                    pctlpts=80 90 
                    pctlname=abtdi80 abtdi90
   ;
run;

data omoa;
   set stats;
   label msd="Mean Squared Deviation"
         mad="Mean Absolute Deviation"
		 mdiff="Mean Difference"
		 P_abtdi80="TDI 80"
		 P_abtdi90="TDI 90"
		 ;
run;

ods select all;

proc print label noobs;
   var mdiff msd mad P_abtdi80 P_abtdi90;
   Title4 "Other Measures of Pairwise Observer Agreement";
run;

ods select none;

%end; *** End pairwise block ***;


%end;  * End NON-REPLICATE Coefficient calculation *;




*** Replication block ***;


%if &numrep ge 2 %then %do; *** Only do this if there are replicates in the data ***;

**** Begin block for calculating CIV ****;

Title4 "Coefficient of Interobserver Variability from Sums of Squares";

data CIV;
    set &lib..&dsn;

    **** set up arrays, y(j,k), yj=y(j,.), u(j);
    array y(&numrater,&numrep) &varlist;
    array yj(&numrater) yj1-yj&numrater;
    array u(&numrater) u1-u&numrater;
 
    ******** initialize ydd ( y..);

    ydd=0;

    ********* calculate v;
    do j=1 to &numrater;
	   yj(j)=0;  ******* this is yj.;
	   do k=1 to &numrep;
	      ydd=ydd+y(j,k);   ***** get ydd.;
		  yj(j)=yj(j)+y(j,k);  **** get yj.;
       end;
	end;

	*** Calculate the means for ydd and yj ***;

	yddbar = ydd/%eval(&numrater*&numrep);

	array yjbar (&numrater) yjbar1-yjbar&numrater ;

	do j=1 to &numrater;
	   yjbar(j)=yj(j)/&numrep ;
	end;

	v=0;

    do j=1 to &numrater;
	   v=v+(yjbar(j)-yddbar)**2;
    end;

    v=v/(&numrater-1);

    *****  calculate u(j) and u. ;
	**** u. is used in calculation of t2 and w2;
    ud=0;

    do j=1 to &numrater;
	   u(j)=0;
	   do k=1 to &numrep;
	      u(j)=u(j)+(y(j,k)-yjbar(j))**2;
	   end;
	   u(j)=u(j)/(&numrep-1);
	   ud=ud+u(j);
	end;

	ud=ud/(&numrater);

    keep v ud u1 - u&numrater;

run;

proc means noprint;
   var v ud u1 - u&numrater;
   output out=CIVa mean=vd udd u1 - u&numrater;  **** all we need are means;
run;

data CIVb;
   set CIVa;
   t2=vd-(udd/&numrep);
   w2=udd;
   civ=t2/(t2+w2);
   psi=1-civ;
   keep vd udd t2 w2 civ psi u1 - u&numrater;
   label vd='V.'
	     udd='U..'
		 t2='Tau Squared'
		 w2='Omega Squared'
		 civ='CIV'
		 psi='PSI'
   ;
run;


*** Bootstrap Intervals for CIV and PSI ***;

** 1. Generate the bootstrap samples, using PRC SURVEYSELECT to get samples wr **;

%if &boot gt 0 %then %do;

ods select none;

proc surveyselect data=&lib..&dsn method=urs n=&numsub out=bootsamp rep=&boot;
run;

data boot1;
   set bootsamp;
   do i=1 to numberhits;
      output;
   end;
   drop i numberhits ;
run;

data CIVboot;
    set boot1;
	by replicate;

    **** set up arrays, y(j,k), yj=y(j,.), u(j);
    array y(&numrater,&numrep) &varlist;
    array yj(&numrater) yj1-yj&numrater;
    array u(&numrater) u1-u&numrater;
 
    ******** initialize ydd ( y..);

    ydd=0;

    ********* calculate v;
    do j=1 to &numrater;
	   yj(j)=0;  ******* this is yj.;
	   do k=1 to &numrep;
	      ydd=ydd+y(j,k);   ***** get ydd.;
		  yj(j)=yj(j)+y(j,k);  **** get yj.;
       end;
	end;

	*** Calculate the means for ydd and yj ***;

	yddbar = ydd/%eval(&numrater*&numrep);

	array yjbar (&numrater) yjbar1-yjbar&numrater ;

	do j=1 to &numrater;
	   yjbar(j)=yj(j)/&numrep ;
	end;

	v=0;

    do j=1 to &numrater;
	   v=v+(yjbar(j)-yddbar)**2;
    end;

    v=v/(&numrater-1);

    *****  calculate u(j) and u. ;
	**** u. is used in calculation of t2 and w2;
    ud=0;

    do j=1 to &numrater;
	   u(j)=0;
	   do k=1 to &numrep;
	      u(j)=u(j)+(y(j,k)-yjbar(j))**2;
	   end;
	   u(j)=u(j)/(&numrep-1);
	   ud=ud+u(j);
	end;

	ud=ud/(&numrater);

    keep v ud u1 - u&numrater replicate;

run;

proc means noprint;
   var v ud u1 - u&numrater;
   by replicate;
   output out=bootCIVa mean=vd udd u1 - u&numrater;  **** all we need are means;
run;

data bootCIVb;
   set bootCIVa;
   t2=vd-(udd/&numrep);
   w2=udd;
   civ=t2/(t2+w2);
   psi=1-civ;
   keep vd udd t2 w2 civ psi u1 - u&numrater;
   label vd='V.'
	     udd='U..'
		 t2='Tau Squared'
		 w2='Omega Squared'
		 civ='CIV'
		 psi='PSI'
   ;
run;

proc univariate data=bootcivb;
   var civ psi;
   output out=bootint pctlpre=c_ p_ 
                      pctlpts=%sysevalf(100*&alphalow) %sysevalf(100*&alphahi) 
                      pctlname=low high ;
run;

data fullciv;
   merge civb bootint;
   label c_low="Lower &conf %"
         c_high="Upper &conf %"
         p_low="Lower &conf %"
         p_high="Upper &conf %"
		 ;
run;

%end; *** Bootstrap block ***;

ods select all;

%if &boot gt 0 %then %do;

proc print label noobs;
   var vd udd t2 w2 civ c_low c_high;
run;

proc print label noobs;
   var vd udd t2 w2 psi p_low p_high;
run;


%end;
%else %do;

proc print label noobs;
   var vd udd t2 w2 civ psi;
run;

%end;

ods select none;

*** posthocs for CIV and PSI ***;

%if &numrater ge 3 %then %do; *** Begin Posthocs block ***;

%do x=1 %to %eval(&numrater - 1);
%do y=%eval(&x+1) %to &numrater;


data civp; *** Create a data set with only two measures, repeat for all pairwise measures ***;
   set &lib..&dsn;
   %do r=1 %to &numrep;

     var1&r=&var&x&r; * Pick first measure *;

   %end;
   %do r=1 %to &numrep;

     var2&r=&var&y&r; * Pick second measure *;

   %end;
   keep id 

   %do r=1 %to &numrep;

   var1&r

   %end; 

   %do r=1 %to &numrep;

   var2&r

   %end; 


   ;; * Keep data set with chosen measures as measure 1 and measure 2 *;
run;

ods select none;

data CIVp2;
    set civp;

    **** set up arrays, y(j,k), yj=y(j,.), u(j);
    array y(2,&numrep) var11 -- var2&numrep;
    array yj(2) yj1 yj2;
    array u(2) u1 u2;
 
    ******** initialize ydd ( y..);

    ydd=0;

    ********* calculate v;
    do j=1 to 2;
	   yj(j)=0;  ******* this is yj.;
	   do k=1 to &numrep;
	      ydd=ydd+y(j,k);   ***** get ydd.;
		  yj(j)=yj(j)+y(j,k);  **** get yj.;
       end;
	end;

	*** Calculate the means for ydd and yj ***;

	yddbar = ydd/%eval(2*&numrep);

	array yjbar (2) yjbar1 yjbar2 ;

	do j=1 to 2;
	   yjbar(j)=yj(j)/&numrep ;
	end;

	v=0;

    do j=1 to 2;
	   v=v+(yjbar(j)-yddbar)**2;
    end;

    v=v/(2-1);

    *****  calculate u(j) and u. ;
	**** u. is used in calculation of t2 and w2;
    ud=0;

    do j=1 to 2;
	   u(j)=0;
	   do k=1 to &numrep;
	      u(j)=u(j)+(y(j,k)-yjbar(j))**2;
	   end;
	   u(j)=u(j)/(&numrep-1);
	   ud=ud+u(j);
	end;

	ud=ud/(2);

    keep v ud u1 u2;

run;

proc means noprint data=CIVp2;
   var v ud u1 u2;
   output out=CIVa mean=vd udd u1 u2;  **** all we need are means;
run;

data CIV&x&y;
   set CIVa;
   rater1 = &x;
   rater2 = &y;
   t2=vd-(udd/&numrep);
   w2=udd;
   civ=t2/(t2+w2);
   psi=1-civ;
   keep rater1 rater2 vd udd t2 w2 civ psi u1 u2;
   label rater1='First Rater'
         rater2='Second Rater'
         vd='V.'
	     udd='U..'
		 t2='Tau Squared'
		 w2='Omega Squared'
		 civ='CIV'
		 psi='PSI'
   ;
run;

*** Bootstrap Intervals for CIV and PSI ***;

** 1. Generate the bootstrap samples, using PRC SURVEYSELECT to get samples wr **;

%if &boot gt 0 %then %do;

ods select none;

proc surveyselect data=civp method=urs n=&numsub out=pbootsamp rep=&boot; *** Bootstrap samples from the posthoc dataset ***;
run;

data boot1;
   set pbootsamp;
   do i=1 to numberhits;
      output;
   end;
   drop i numberhits ;
run;

data CIVboot;
    set boot1;
	by replicate;

    **** set up arrays, y(j,k), yj=y(j,.), u(j);
    array y(2,&numrep) var11 -- var2&numrep;
    array yj(2) yj1 yj2;
    array u(2) u1 u2;
 
    ******** initialize ydd ( y..);

    ydd=0;

    ********* calculate v;
    do j=1 to 2;
	   yj(j)=0;  ******* this is yj.;
	   do k=1 to &numrep;
	      ydd=ydd+y(j,k);   ***** get ydd.;
		  yj(j)=yj(j)+y(j,k);  **** get yj.;
       end;
	end;

	*** Calculate the means for ydd and yj ***;

	yddbar = ydd/%eval(2*&numrep);

	array yjbar (2) yjbar1 yjbar2 ;

	do j=1 to 2;
	   yjbar(j)=yj(j)/&numrep ;
	end;

	v=0;

    do j=1 to 2;
	   v=v+(yjbar(j)-yddbar)**2;
    end;

    v=v;

    *****  calculate u(j) and u. ;
	**** u. is used in calculation of t2 and w2;
    ud=0;

    do j=1 to 2;
	   u(j)=0;
	   do k=1 to 2;
	      u(j)=u(j)+(y(j,k)-yjbar(j))**2;
	   end;
	   u(j)=u(j);
	   ud=ud+u(j);
	end;

	ud=ud/(2);

    keep v ud u1 u2 replicate;

run;

proc means noprint;
   var v ud u1 u2;
   by replicate;
   output out=bootCIVa mean=vd udd u1 u2;  **** all we need are means;
run;

data bootCIVb;
   set bootCIVa;
   t2=vd-(udd/2);
   w2=udd;
   civ=t2/(t2+w2);
   psi=1-civ;
   keep vd udd t2 w2 civ psi u1 u2;
   label vd='V.'
	     udd='U..'
		 t2='Tau Squared'
		 w2='Omega Squared'
		 civ='CIV'
		 psi='PSI'
   ;
run;

proc univariate data=bootcivb;
   var civ psi;
   output out=bootint pctlpre=c_ p_ 
                      pctlpts=%sysevalf(100*&alphalow) %sysevalf(100*&alphahi) 
                      pctlname=low high ;
run;

data fullciv&x&y;
   merge civ&x&y bootint;
   rater1=&x;
   rater2=&y;
   label c_low="Lower &conf %"
         c_high="Upper &conf %"
         p_low="Lower &conf %"
         p_high="Upper &conf %"
		 ;
run;

%end; *** Bootstrap block ***;





%end;
%end;

data civph;
   set

   %if &boot eq 0 %then %do;

   %do x=1 %to %eval(&numrater - 1);
       %do y=%eval(&x+1) %to &numrater;
	       civ&x&y
	   %end;
   %end;

   ;;

   run;

ods select all;

Title5 "Pairwise By Rater";

proc print noobs data=civph label;
   var rater1 rater2 civ psi;
run;


   %end; %else %do;

   %do x=1 %to %eval(&numrater - 1);
       %do y=%eval(&x+1) %to &numrater;
	       fullciv&x&y
	   %end;
   %end;

   ;;

   run;

ods select all;

Title5 "Pairwise By Rater";

proc print noobs data=civph label;
   var rater1 rater2 civ c_low c_high psi p_low p_high;
run;



   %end;

run;

ods select none;

%end;

**** Moments and CCCs, Method 1 ****;

Title4 "Moments and CCCs, Method 1";
title5 "Based on the First Observation of each Observer";

data mccc1;

   set &lib..&dsn;

   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      frst(j)=y(j,1);
   end;

   keep f1 - f&numrater;

run;

proc means noprint data=mccc1;
   var f1 - f&numrater;
   output out=mccc1a mean=m1 - m&numrater
                     var =v1 - v&numrater
   ;
run;

data mccc1b;
   merge mccc1a CIVb;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater;

      w(j) = u(j);
	  d(j) = sigma(j) - w(j);

   end;

run;

data mccc1c;
   set mccc1b;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater ;

      rater=j;
      muhat=m(j);
	  deltahatsquared=d(j);
	  whatsquared=w(j);
	  sigmahatsquared=sigma(j);
	  sigmahat=sqrt(sigma(j));

	  output;

   end;

   label rater="Rater"
         muhat="Mu"
		 deltahatsquared="Delta Squared"
		 whatsquared="Omega Squared"
		 sigmahatsquared="Sigma Squared"
		 sigmahat="Sigma"
   ;

run;

ods select none; *** Want to keep CORR matrix, but not output to screen ***;

*** Suppress printing method 1 results until further notice ***;

proc print noobs label;
   var rater muhat deltahatsquared whatsquared sigmahatsquared sigmahat;
run;
 
proc corr data=mccc1 ;
   var f1 - f&numrater;
   ods output PearsonCorr=corrs;
run;

%do d=1 %to &numrater;

data r&d;
   set corrs;
   if _n_=&d;
   %do j=1 %to &numrater;
   r&d&j=f&j;
   %end;
run;

%end;

data r1line;
   merge
   %do d=1 %to &numrater;
   r&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       r&i&j
       %end;
   %end;
   ;;

run;

data mccc1d;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;

	     rater1 = j;
		 rater2 = jprime;

		 rhohat=rho(j,jprime);

		 rhohatmu = rho(j,jprime) * sqrt( sigma(j)*sigma(jprime) / (d(j)*d(jprime)) ) ;

		 output;

	  end;
   end;


   label rater1="Rater 1"
         rater2="Rater 2"
		 rhohat="Pearson Correlation of First Measures"
		 rhohatmu="Rho-mu"
   ;
   
run;

data mccc1f;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;
   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   denompij=0;

   do j=1 to %eval(&numrater-1);
      do jprime = j+1 to &numrater;
	     denompij=denompij + (rho(j,jprime) * sqrt(sigma(j) * sigma(jprime)) / 
                             (&numrater * (&numrater-1)/2));
	  end;
   end;

   do j=1 to &numrater;
      rater=j;
	  mu_j=m(j);
	  delta_sq_j = d(j);
	  omega_sq_j = w(j);
	  sigma_sq_j = sigma(j);
      rhohatij= d(j)/sigma(j) ;
	  pij = sigma(j) / denompij;
	  checker=pij*(1-rhohatij);
	  output;
   end;

   label mu_j="Mu"
         delta_sq_j="Delta- Squared"
         omega_sq_j="Omega- Squared"
         sigma_sq_j="Sigma- Squared"
         rhohatij="Rho^I"
		 pij="Pi"
		 checker="Product for checking rho-hats"
		 ;

run;



**** Overall Measures of Agreement, Method 1 ****;
**** Calculated using first replicate        ****;

**** Use code from non-rep (above) and modify to handle first obs from data ****;

data cccr;
   set &lib..&dsn;
   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      frst(j)=y(j,1);
   end;

   keep f1 - f&numrater;

run;

ods select none;
proc corr cov;
   var f1-f&numrater;
   ods output cov=covmat;
run;

%do d=1 %to &numrater;

data a&d;
   set covmat;
   if _n_=&d;
   %do j=1 %to &numrater;
   &var&d&j=f&j;
   %end;
run;

%end;

data oneline;
   merge
   %do d=1 %to &numrater;
   a&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       &var&i&j
       %end;
   %end;
   ;; 
run;

data pear;
   set oneline;

   %do i=1 %to &numrater;
       %do j=1 %to &numrater;

	   r&i&j = &var&i&j / sqrt(&var&i&i * &var&j&j);

	   %end;
   %end;

run;

%let a11 = 11;

data corrs;
   set pear;
   array r(&numrater,&numrater) r11 -- r&numrater&numrater;
   array x(&numrater,&numrater) &var&a11 -- &var&numrater&numrater;
   do i=1 to &numrater;
      do j=1 to &numrater;
	     pearson=r(i,j);
		 cov=x(i,j);
		 var1=x(i,i);
		 var2=x(j,j);
		 if i lt j then output;
	  end;
   end;
   keep i j pearson cov var1 var2;
   label pearson="Pearson Correlation"
         cov="Covariance"
         var1="Variance 1"
		 var2="Variance 2"
		 i="Measure 1"
		 j="Measure 2"
		 ;
run;

*** Calculating the OCCC ***;

* Need Covariances from Pearson *;

proc means sum noprint data=corrs;
   var cov;
   output out=covs sum=sumcov;
run;

proc sort data=corrs;
   by i j;
run;

data v1;
   set corrs;
   by i j;
   if first.i; *** Keeps first j-1 variances ***;
   measure=i;
   var=var1;
   keep measure var;
run;
   
data v2;
   set corrs;
   by i j;
   if last.j; *** Keeps last j-1 variances ***;
   measure=j;
   var=var2;
   keep measure var;
run;

data v;
   set v1 v2;
run;

proc sort;
   by measure;
run;

data vx;
   set v;
   by measure;
   if first.measure;
run;

proc means sum noprint data=vx;
   var var;
   output out=vars sum=sumvar;
run;

proc means mean data=cccr noprint;
   var f1-f&numrater;
   output out=means mean=;
run;

data meandiffs;
   set means;
   array means (&numrater) f1-f&numrater;
   xend = &numrater - 1;
   meansq=0; * to start *;
   do i = 1 to xend;
      do j = i+1 to &numrater;
	     meansq + (means(i) - means(j))**2;
	  end;
   end;
   drop i j xend;
run;

proc means mean data=corrs noprint;
   var pearson;
   output out=meancorr mean=corrbar;
run;

data occc;
   merge covs vars meandiffs meancorr;
   nab = 2*sumcov / (((&numrater - 1)*sumvar) +  meansq);
   label nab = "Coefficient of Absolute Agreement"
		 ;
run;


data mccc1e;
   merge mccc1b r1line occc;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;
   num=0;

   rhohatmuc=0;
   muss=(
   %do k = 1 %to %eval(&numrater-1);
   m&k +
   %end;
   m&numrater)/&numrater;
   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;
         num=num+2*(sqrt((sigma(j)*sigma(jprime))*rho(j,jprime)));
	  end;
   end;
   denom=((&numrater*(
		 %do i=1 %to %eval(&numrater-1);
         ((m&i - muss)**2) +
		 %end;
		 (m&numrater - muss)**2) +
		 (&numrater-1)*(
		 %do i=1 %to %eval(&numrater-1);
         d&i +
		 %end;
         d&numrater) ));
   rhohatmuc=num/denom;
   
   label rhohatmuc="Rho-c (mu)"
         nab="Rho-c (y)"
		 ;    

run;

*** Check sum ***;

data check;
   set mccc1f;
   cs+(pij * (1 - rhohatIJ)/&numrater);
   if _n_ = &numrater;
   keep cs;
run;
   
data check2;
   merge mccc1e check;
   left=1/nab;
   right=(1/rhohatmuc) +cs;
   keep left right;
run;

**** Now we have a data set with the Ordinary CCCs ****;

*ods select all;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print noobs data=mccc1f label;
   var rater mu_j delta_sq_j omega_sq_j sigma_sq_j pij rhohatij checker;
run;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print data=mccc1e label noobs;
   var nab rhohatmuc;
run;

proc print data=check2 noobs;
   var left right;
   title6 "Checking Rho-c(y) and Rho-c(mu)";
run;


ods select none;


**** END METHOD 1 BLOCK ****;






**** Moments and CCCs, Method 2 ****;

Title4 "Moments and CCCs, Method 2";
title5 "Based on a Randomly Selected Observation from each Observer";

data mccc1;

   set &lib..&dsn;

   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      rand=rantbl(&seed,
	  %do k=1 %to (&numrep-1);
	  1/&numrep,
	  %end;
      1/&numrep);
      frst(j)=y(j,rand);
   end;

   *keep f1 - f&numrater;

run;

proc means noprint;
   var f1 - f&numrater;
   output out=mccc1a mean=m1 - m&numrater
                     var =v1 - v&numrater
   ;
run;

data mccc1b;
   merge mccc1a CIVb;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater;

      w(j) = u(j);
	  d(j) = sigma(j) - w(j);

   end;

run;

data mccc1c;
   set mccc1b;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater ;

      rater=j;
      muhat=m(j);
	  deltahatsquared=d(j);
	  whatsquared=w(j);
	  sigmahatsquared=sigma(j);
	  sigmahat=sqrt(sigma(j));

	  output;

   end;

   label rater="Rater"
         muhat="Mu"
		 deltahatsquared="Delta Squared"
		 whatsquared="Omega Squared"
		 sigmahatsquared="Sigma Squared"
		 sigmahat="Sigma"
   ;

run;

ods select none; *** suppress printing for method 2 until further notice ***;

proc print noobs label;
   var rater muhat deltahatsquared whatsquared sigmahatsquared sigmahat;

run;
 
ods select none; *** Want to keep CORR matrix, but not output to screen ***;
proc corr data=mccc1 ;
   var f1 - f&numrater;
   ods output PearsonCorr=corrs;
run;

%do d=1 %to &numrater;

data r&d;
   set corrs;
   if _n_=&d;
   %do j=1 %to &numrater;
   r&d&j=f&j;
   %end;
run;

%end;

data r1line;
   merge
   %do d=1 %to &numrater;
   r&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       r&i&j
       %end;
   %end;
   ;;

run;

data mccc1d;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;

	     rater1 = j;
		 rater2 = jprime;

		 rhohat=rho(j,jprime);

		 rhohatmu = rho(j,jprime)*sqrt(sigma(j)*sigma(jprime)/(d(j)*d(jprime))) ;

		 output;

	  end;
   end;

   label rater1="Rater 1"
         rater2="Rater 2"
		 rhohat="Pearson Correlation of First Measures"
		 rhohatmu="Rho-mu"
   ;

run;

data mccc1f;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;
   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   denompij=0;

   do j=1 to %eval(&numrater-1);
      do jprime = j+1 to &numrater;
	     denompij=denompij + (rho(j,jprime) * sqrt(sigma(j) * sigma(jprime)) / 
                             (&numrater * (&numrater-1)/2));
	  end;
   end;

   do j=1 to &numrater;
      rater=j;
	  mu_j=m(j);
	  delta_sq_j = d(j);
	  omega_sq_j = w(j);
	  sigma_sq_j = sigma(j);
      rhohatij= d(j)/sigma(j) ;
	  pij = sigma(j) / denompij;
	  checker=pij*(1-rhohatij);
	  output;
   end;

   label mu_j="Mu"
         delta_sq_j="Delta- Squared"
         omega_sq_j="Omega- Squared"
         sigma_sq_j="Sigma- Squared"
         rhohatij="Rho^I"
		 pij="Pi"
		 checker="Product for checking rho-hats"
		 ;

run;


data cccr;
   set &lib..&dsn;
   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      frst(j)=y(j,1);
   end;

   keep f1 - f&numrater;

run;

ods select none;
proc corr cov;
   var f1-f&numrater;
   ods output cov=covmat;
run;

%do d=1 %to &numrater;

data a&d;
   set covmat;
   if _n_=&d;
   %do j=1 %to &numrater;
   &var&d&j=f&j;
   %end;
run;

%end;

data oneline;
   merge
   %do d=1 %to &numrater;
   a&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       &var&i&j
       %end;
   %end;
   ;; 
run;

data pear;
   set oneline;

   %do i=1 %to &numrater;
       %do j=1 %to &numrater;

	   r&i&j = &var&i&j / sqrt(&var&i&i * &var&j&j);

	   %end;
   %end;

run;

%let a11 = 11;

data corrs;
   set pear;
   array r(&numrater,&numrater) r11 -- r&numrater&numrater;
   array x(&numrater,&numrater) &var&a11 -- &var&numrater&numrater;
   do i=1 to &numrater;
      do j=1 to &numrater;
	     pearson=r(i,j);
		 cov=x(i,j);
		 var1=x(i,i);
		 var2=x(j,j);
		 if i lt j then output;
	  end;
   end;
   keep i j pearson cov var1 var2;
   label pearson="Pearson Correlation"
         cov="Covariance"
         var1="Variance 1"
		 var2="Variance 2"
		 i="Measure 1"
		 j="Measure 2"
		 ;
run;

*** Calculating the OCCC ***;

* Need Covariances from Pearson *;

proc means sum noprint data=corrs;
   var cov;
   output out=covs sum=sumcov;
run;

proc sort data=corrs;
   by i j;
run;

data v1;
   set corrs;
   by i j;
   if first.i; *** Keeps first j-1 variances ***;
   measure=i;
   var=var1;
   keep measure var;
run;
   
data v2;
   set corrs;
   by i j;
   if last.j; *** Keeps last j-1 variances ***;
   measure=j;
   var=var2;
   keep measure var;
run;

data v;
   set v1 v2;
run;

proc sort;
   by measure;
run;

data vx;
   set v;
   by measure;
   if first.measure;
run;

proc means sum noprint data=vx;
   var var;
   output out=vars sum=sumvar;
run;

proc means mean data=cccr noprint;
   var f1-f&numrater;
   output out=means mean=;
run;

data meandiffs;
   set means;
   array means (&numrater) f1-f&numrater;
   xend = &numrater - 1;
   meansq=0; * to start *;
   do i = 1 to xend;
      do j = i+1 to &numrater;
	     meansq + (means(i) - means(j))**2;
	  end;
   end;
   drop i j xend;
run;

proc means mean data=corrs noprint;
   var pearson;
   output out=meancorr mean=corrbar;
run;

data occc;
   merge covs vars meandiffs meancorr;
   nab = 2*sumcov / (((&numrater - 1)*sumvar) + meansq);
   label nab = "Coefficient of Absolute Agreement"
		 ;
run;


data mccc1e;
   merge mccc1b r1line occc;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;
   num=0;

   rhohatmuc=0;
   muss=(
   %do k = 1 %to %eval(&numrater-1);
   m&k +
   %end;
   m&numrater)/&numrater;
   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;
         num=num+2*(sqrt((sigma(j)*sigma(jprime))*rho(j,jprime)));
	  end;
   end;
   denom=((&numrater*(
		 %do i=1 %to %eval(&numrater-1);
         ((m&i - muss)**2) +
		 %end;
		 (m&numrater - muss)**2) +
		 (&numrater-1)*(
		 %do i=1 %to %eval(&numrater-1);
         d&i +
		 %end;
         d&numrater) ));
   rhohatmuc=num/denom;
   
   label rhohatmuc="Rho-c (mu)"
         nab="Rho-c (y)"
		 ;    

run;

data check;
   set mccc1f;
   cs+pij * (1 - rhohatIJ)/&numrater;
   if _n_ = &numrater;
run;
   
data check2;
   merge mccc1e check;
   left=1/nab;
   right=1/rhohatmuc+cs;
   keep left right;
run;

ods select none; *** suppress printing for method 2 until further notice ***;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print noobs data=mccc1f label;
   var rater mu_j delta_sq_j omega_sq_j sigma_sq_j pij rhohatij checker;
run;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print data=mccc1e label noobs;
   var nab rhohatmuc;
run;

proc print data=check2 noobs;
   var left right;
   title6 "Checking Rho-c(y) and Rho-c(mu)";
run;

ods select none;


*** Using proc mixed ***;


data mixed;
   set &lib..&dsn; *** This is the GLM data set ***;

   array mix (&numrater,&numrep) &varlist  ; 

   do rater=1 to &numrater;
      do rep=1 to &numrep; 

         &var=mix(rater, rep); 
         output;	     

	  end;
   end;
   keep &var rater id;
run;

proc mixed data=mixed;
   class id rater ;
   model &var=rater / noint solution;
   random rater / g subject=id type=un v;
   repeated / r group=rater;
   title;
   ods output g=gmatrix covparms=ests solutionf=means;
run;

data resids;
   set ests;
   if covparm="Residual"; *** Keeps sigma-squares only ***;

   %do k=1 %to &numrater;
   if group="rater &k" then do;
      effect="rater"; rater=&k; sigsq=estimate;
   end; 
   %end;
   drop estimate;
run;
proc print;
run;

data mus;
   set means;
   mu=estimate;
   keep effect rater mu;
run;

data all;
   merge gmatrix resids mus;
   by effect rater;

   array g (&numrater) 
   %do k=1 %to &numrater;
   col&k
   %end;
   ;;

   do i=1 to &numrater; *** Iterate over raters ***;

      if rater=i then deltasq=g(i); *** Identify the deltas ***;

   end; 

   rho_i = deltasq/(deltasq+sigsq); *** Calculate rho^I for each rater ***;

   drop covparm subject;

run;
ods select all;
proc print;
   var rater deltasq sigsq rho_i;
   title5 "Agreement coefficients from the Mixed Model parameters";
run;
ods select none;

%do k=1 %to &numrater;

data rater&k;
   set all;
   if rater=&k;
   ssq&k=sigsq;
   %do kp=1 %to &numrater;
       rho&k&kp=col&kp; 
   %end;
   mu&k=mu; rho&k=rho_i;
   keep ssq&k -- rho&k;
run;

%end;

data oneline;
   merge 
   %do k=1 %to &numrater; 
         rater&k
   %end;
   ;;

   *** Now, all of the data are on one line -- we can use arrays! ***;

   array rho (&numrater,&numrater) 
   %do k=1 %to &numrater;
       %do kp=1 %to &numrater;
	       rho&k&kp
	   %end;
   %end;
   ;;

   array dsq (&numrater) 
   %do k=1 %to &numrater;
	       rho&k&k
   %end;
   ;;

   array ssq (&numrater) 
   %do k=1 %to &numrater;
	       ssq&k
   %end;
   ;;

   array mu (&numrater)
   %do k=1 %to &numrater;
	       mu&k
   %end;
   ;;


   do j=1 to &numrater;

      sumdsq+dsq(j);
	  sumsq+ssq(j);

   end;

   do j=1 to &numrater-1;
      do jp=j+1 to &numrater;

	     num+rho(j,jp);
		 den+((mu(j) - mu(jp))**2);

	  end;
   end;

   rho_c_mu = (2*num)/((&numrater-1)*sumdsq+den);
   rho_c_y = (2*num)/((&numrater-1)*sumdsq+(&numrater-1)*sumsq+den);

run;

ods select all;
proc print noobs;
   var rho_c_mu rho_c_y;

run;
ods select none;


*** posthocs for Rho-mu and Rho-c ***;

%if &numrater ge 3 %then %do; *** Begin Posthocs block ***;

%do x=1 %to %eval(&numrater - 1);
%do y=%eval(&x+1) %to &numrater;

*** Using proc mixed ***;


data m; *** Create a data set with only two measures, repeat for all pairwise measures ***;
   set &lib..&dsn;
   %do r=1 %to &numrep;

     var1&r=&var&x&r; * Pick first measure *;

   %end;
   %do r=1 %to &numrep;

     var2&r=&var&y&r; * Pick second measure *;

   %end;
   keep id 

   %do r=1 %to &numrep;

   var1&r

   %end; 

   %do r=1 %to &numrep;

   var2&r

   %end; 


   ;; * Keep data set with chosen measures as measure 1 and measure 2 *;
run;





data mixed;
   set m; *** This is the GLM data set ***;

   array mix (2,&numrep) 
   %do r=1 %to &numrep;

   var1&r

   %end; 

   %do r=1 %to &numrep;

   var2&r

   %end; 
;; 

   do rater=1 to 2;
      do rep=1 to &numrep; 

         &var=mix(rater, rep); 
         output;	     

	  end;
   end;
   keep &var rater id;
run;

ods select none;  *** Checking to ensure that the data look right ***;
proc print data=mixed (obs=40);
title "Raters &x and &y";
run;
ods select none;

proc mixed data=mixed;
   class id rater ;
   model &var=rater / noint solution;
   random rater / g subject=id type=un v;
   repeated / r group=rater;
   title;
   ods output g=gmatrix covparms=ests solutionf=means;
run;

data resids;
   set ests;
   if covparm="Residual"; *** Keeps sigma-squares only ***;

   %do k=1 %to 2;
   if group="rater &k" then do;
      effect="rater"; rater=&k; sigsq=estimate;
   end; 
   %end;
   drop estimate;
run;
proc print;
run;

data mus;
   set means;
   mu=estimate;
   keep effect rater mu;
run;

data all;
   merge gmatrix resids mus;
   by effect rater;

   array g (2) 
   %do k=1 %to 2;
   col&k
   %end;
   ;;

   do i=1 to 2; *** Iterate over raters ***;

      if rater=i then deltasq=g(i); *** Identify the deltas ***;

   end; 

   rho_i = deltasq/(deltasq+sigsq); *** Calculate rho^I for each rater ***;

   drop covparm subject;

run;
ods select none;
proc print;
   var rater deltasq sigsq rho_i;
   title5 "Agreement coefficients from the Mixed Model parameters";
run;
ods select none;

%do k=1 %to 2;

data rater&k;
   set all;
   if rater=&k;
   ssq&k=sigsq;
   %do kp=1 %to 2;
       rho&k&kp=col&kp; 
   %end;
   mu&k=mu; rho&k=rho_i;
   keep ssq&k -- rho&k;
run;

%end;

data oneline&x&y;
   merge 
   %do k=1 %to 2; 
         rater&k
   %end;
   ;;

   *** Now, all of the data are on one line -- we can use arrays! ***;

   array rho (2,2) 
   %do k=1 %to 2;
       %do kp=1 %to 2;
	       rho&k&kp
	   %end;
   %end;
   ;;

   array dsq (2) 
   %do k=1 %to 2;
	       rho&k&k
   %end;
   ;;

   array ssq (2) 
   %do k=1 %to 2;
	       ssq&k
   %end;
   ;;

   array mu (2)
   %do k=1 %to 2;
	       mu&k
   %end;
   ;;


   do j=1 to 2;

      sumdsq+dsq(j);
	  sumsq+ssq(j);

   end;

   do j=1 to 2-1;
      do jp=j+1 to 2;

	     num+rho(j,jp);
		 den+((mu(j) - mu(jp))**2);

	  end;
   end;

   rho_c_mu = (2*num)/(sumdsq+den);
   rho_c_y = (2*num)/(sumdsq+sumsq+den);
   rater1 = &x;
   rater2 = &y;

run;

%end; *** Observer 2 loop ***;
%end; *** Observer 1 loop ***;

data rhos;
set
%do x=1 %to %eval(&numrater - 1);
%do y=%eval(&x+1) %to &numrater;
oneline&x&y
%end;
%end;

;;
run;

ods select all;
proc print data=rhos noobs;
   var rater1 rater2 rho_c_mu rho_c_y;
run;
ods select none;



%end; *** Pairwise loop ***;



**** END METHOD 2 BLOCK ****;





**** Overall Measures of Agreement from ANOVA model ****;

title4 "Measures of Agreement based on 2-way Mixed Model with Interaction";

data anovaprep;
   set &lib..&dsn;

   %do j=1 %to &numrater; *** repeat over raters ***;

       %do k=1 %to &numrep; *** repeat over replicates ***;

	       y = &var&j&k; rater=&j; rep=&k; subject=id; output; *** create a dataset with each measure on its
		                                                          *** own line. ;

	   %end; *** replicate loop ***;

   %end; *** rater loop ***;

   keep rater subject y rep;

run;
 
      ************** anova method starts here;
ods select none;

proc glm data=anovaprep;
   class subject rater;
   model y=rater subject subject*rater  ;
   random subject / test;
   ods output modelanova=modelanova;
run;

data anovaout;
   set modelanova;
   if Hypothesistype=3;
run;

data sr;
   set anovaout;
   if source='rater' ;
   msr=ms;
   keep msr;
run;

data ss;
    set anovaout;
	if source='subject' ;
	mss=ms;
	keep mss;
run;

data ssr;
    set anovaout;
	if source='subject*rater' ;
	mssr=ms;
	mse=ms/fvalue;
	keep mssr mse;
run;

data combo;
   merge sr ss ssr; 
   ss=(mss-mse)/(&numrep*&numrater);
   sr=(msr-mssr)/(&numrep*&numsub);
   ssr=(mssr-mse)/&numrep;
   se=mse;
   civ=(sr+ssr)/(sr+ssr+se);
   rhoy=(ss-(ssr/(&numrater-1)))/(ss+sr+ssr+se);
   rhou=(ss-(ssr/(&numrater-1)))/(ss+sr+ssr);
   keep ss sr ssr se civ rhoy rhou;
   label civ='CIV'
         rhoy='Rho-c (y) From Variance Components'
         rhou='Rho-c (mu) From Variance Components'
		 ss="Subject Component"
		 sr="Rater Component"
		 ssr="Subject- Rater Interaction Component"
		 se="Error Component"
   ;
run;

ods select all;
proc print noobs label;
   var ss sr ssr se civ rhoy rhou;
run;

%end;  *** End block for replicated measures ***;

ods select all;

%mend ICCs;